Quantum Disorder and Quantum Chaos in Andreev Billiards 
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We investigate the crossover from the semiclassical to the quantum description of electron energy 
states in a chaotic metal grain connected to a superconductor. We consider the influence of scattering 
off point impurities (quantum disorder) and of quantum diffraction (quantum chaos) on the electron 
density of states. We show that both the quantum disorder and the quantum chaos open a gap near 
the Fermi energy. The size of the gap is determined by the mean free time in disordered systems 
and by the Ehrenfest time in clean chaotic systems. Particularly, if both times become infinitely 
large, the density of states is gapless, and if either of these times becomes shorter than the electron 
escape time, the density of states is described by random matrix theory. Using the Usadel equation, 
we also study the density of states in a grain connected to a superconductor by a diffusive contact. 
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I. INTRODUCTION 

Condensed matter physics usually deals with quantum 
disorder. Quantum disordered systems contain random 
impurities with the characteristic size smaller than the 
electron Fermi wavelength Af- In these systems, elec- 
trons travel between scattering off impurities along clas- 
sical trajectories, but the scattering is described by quan- 
tum mechanics. Averaging over impurity configurations 
leads to quantum effects in disordered systems. 1 

Recently, new quantum systems - quantum dots - were 
created. 2 In quantum dots the size of irregularities is 
much larger than the electron wavelength, therefore elec- 
tron motion may be determined by a semiclassical theory. 

In the semiclassical limit certain thermodynamic and 
kinetic quantities can be calculated by expansion in pow- 
ers of the Plank's constant Ti. 3 This expansion corre- 
sponds to a series in powers of Af /L, where L is the char- 
acteristic length scale of the system and Af is the electron 
Fermi wave length. For other quantities such semiclassi- 
cal description is justified only for times shorter than the 
Ehrenfest time tE • At longer times a crossover from semi- 
classical to quantum description of these quantities takes 
place. [Ehrenfest 4 demonstrated that the semiclassical 
description of motion of a quantum particle is justified 
for short times.] 

What is the Ehrenfest time in classically chaotic sys- 
tems, where the classical motion is governed by ex- 
ponential divergence of trajectories? It was shown in 
Refs. [5,6,7], that the Ehrenfest time logarithmically de- 
pends on h. Indeed, in semiclassical theory, a particle is 
represented by a wave packet - a superposition of par- 
ticle wave functions. The size of the wave packet can- 
not be smaller than the electron wavelength Af and can- 
not be larger than the characteristic scale L of varia- 
tions of the potential energy. In chaotic systems the 
size of any wave packet exponentially increases in time. 
Thus, the wave packet increases from its initial size Af 
to size d(t) — Af exp(Ai) as time t increases and ex- 
pands to scale of the system L after the Ehrenfest time 
<e = A -1 hxL/Xp, where A is a typical value of the Lya- 
punov exponent. 



We refer to a particle motion at time longer than the 
Ehrenfest time as a regime of quantum chaos. Quantum 
chaos affects different phenomena of condensed matter 
physics. In a metal with random long-range potential 
the weak localization correction to the a. c. -conductance 
at frequency u> is absent at high frequencies, ivt-E ^> 1. 
Nonetheless, at lower frequencies, ajt^ <§C 1, the weak lo- 
calization correction to the conductance reaches its uni- 
versal quantum limit. 

The statistics of electron levels is correctly described 
by the Gutzwiller trace formula at large energies, et~E ^> 
1. At small energies, et^ <C 1, energy levels obey 
the Wigner-Dyson statistics. At intermediate energies, 
st-E ~ 1, the weak localization effects are important for 
understanding of the level statistics. 8 

The shot noise through a chaotic quantum dot is ab- 
sent, if the electron escape time r osc from the dot is 
shorter than the Ehrenfest time, t csc <C £e- In the op- 
posite limit of long t csc ^> t^, the shot noise reaches its 
universal value. 9,10 

The semiclassical method in theory of superconduc- 
tivity was developed by de Gennes and Tinkham 11 and 
Shapoval 12 for the description of the diffusive electron 
scattering at surfaces. The applications 11,12,13 of the 
semiclassical method gave a basis for a belief 14,15 that 
in order to obtain electron properties of a superconduc- 
tor, first, electron states in the superconductor can be 
found as solutions of semiclassical equations and, then, 
these solutions can be averaged over different classical 
trajectories. However, a more detailed analysis demon- 
strated that in certain problems [e.g. in calculations of 
the magnetic penetration depth as a function of an ap- 
plied magnetic field], the semiclassical method does not 
work, and the averaging over classical trajectories gives a 
different result than that obtained by the standard tech- 
nique of averaging over impurities. 16 

According to ref. [5], even for a system with smooth 
long-range disorder, L ^> Ap, the semiclassical method is 
applicable only at time smaller than the Ehrenfest time 
tE- Paper [5] explains the discrepancy between semiclas- 
sical and quantum methods in theory of superconductiv- 
ity: within the semiclassical approach electrons and holes 
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FIG. 1: A sketch of an Andreev billiard: a metal grain (N), 
connected to a superconductor (S). The size of the grain is L, 
and the size of the contact is b. 



move in opposite directions along identical trajectories 
inside a superconductor. Because quantum diffraction 
switches electrons and holes between different trajecto- 
ries, such semiclassical description breaks down at times 
larger than the Ehrcnfest time tg. The Ehrcnfest time t~E 
and the Lyapunov exponent A were calculated for elec- 
tron systems with smooth impurities of the size L 3> Ap 
in ref. [5] . 

In this paper we study the density of states in an An- 
dreev billiard. An Andreev billiard is a small chaotic 
metal grain or semiconductor quantum dot 2 connected 
to a superconductor through a (SN) contact. Electron 
reflection at the contact is described by the Andreev re- 
flection mechanism, 17 when an electron is reflected as 
a hole, moving in the opposite direction. The den- 
sity of states in Andreev billiards was studied both by 
semiclassical 18,19,20 and quantum methods. 20 The semi- 
classical method gives an exponentially small but non- 
zero tail of the density of states for any finite energy, see 
refs. [18,19]. On the contrary, quantum calculations 20 
within random matrix theory predict a gap in the energy 
spectrum. 21 [Both curves of the density of states are pre- 
sented below in Sec. Ill B 2.] 

The difference between semiclassical and random ma- 
trix results for the density of states has attracted theo- 
retical interest to this problem. 22,23 ' 24 Based on the ap- 
proach of refs. [7,25], the authors of ref. [22] gave a qual- 
itative explanation of the difference between the semi- 
classical and the random matrix methods. They argued 
that the crossover happens at energy of the order of the 
inverse Ehrenfest time, 1/te- Papers [23,24] presented 
quantitative results for the energy gap and we discuss 
these results in Sec. V. 

We show that in both disordered and chaotic Andreev 
billiards the density of states has a gap and a square root 
singularity above the gap. We study a disordered An- 
dreev billiard, containing isotropic point scatters, charac- 
terized by arbitrary mean free path Iq. We investigate the 
crossover from the semiclassical 18,19 to quantum regime 20 
as parameter lo/l esc decreases (here Z osc = t osc uf is the 
average length of trajectories in the dot and i>f is the 
Fermi velocity). We further demonstrate that the ran- 
dom matrix result of ref. [20] for reflectionless contact 
is applicable only at intermediate impurity concentra- 
tion, when the mean free path Iq is smaller than a typ- 



ical length of trajectories in the dot, i>ft csc , and longer 
than size b of the SN contact, b <C Iq <C vft csc . As 
the mean free path becomes shorter than the SN contact 
size, l -C b, we find the density of states, described by a 
solution of the Usadel equation. 26 

For the clean chaotic billiards, we apply general meth- 
ods developed in refs. [7,8,9,25] to study the density of 
states. We investigate the crossover from the semiclassi- 
cal limit 18 ' 19 , t^/T csc — > oo, to a random matrix limit, 20 
^e/tosc — * 0, and calculate the energy gap E g as a func- 
tion of fE/Vesc- We find that the energy dependence of 
the density of states has an oscillating component if the 
Ehrenfest time is finite, see also ref. [23] . 

We demonstrate that the density of states has different 
shapes for systems with different values of the strength 
of disorder or chaos. In all cases, the density of states 
has a finite energy gap. The corresponding curves of 
the density of states are different and cannot be trans- 
formed into each other by change of energy scales. Par- 
ticularly, we consider the high energy asymptote of the 
density of states, et csc » 1, and show that the crossover 
between the semiclassical and quantum regimes can be 
understood as a suppression of repetitive trajectories. 

This paper is organized as follows. In the next section 
we present basic parameters of the model and introduce 
the Eilenberger equation. In section III we consider sys- 
tems with quantum disorder. In section IV we study the 
effect of quantum chaos on the density of states. Section 
V contains discussion and conclusions. 



II. MODEL 

We consider a small (i-dimensional metal grain con- 
nected to a superconductor, see Fig. 1. The characteris- 
tic size of the grain is L and its volume is V g oc L d , with 
d being dimension of a system (in two dimensional sys- 
tem V g oc I? stands for the dot's area). We assume that 
the size b of the contact between the metal grain and the 
superconductor is small, so that L > 6, and denote the 
(d — 1) dimension area of the contact by S c oc & d_1 . We 
also include a random clastic potential into our model, 
which determines the effect of impurities or irregulari- 
ties inside the grain and at its boundaries. The length 
scale between consequent scatterings off this potential is 
referred to as the mean free path = vfTq and To is 
the mean free time. According to the relation between 
mean free path Iq and the size of the system L we can 
distinguish ballistic, l >• L, or diffusive, Iq <C L, limits. 

An important energy parameter for a metal grain is 
the Thouless energy, which is inversely proportional to 
the ergodic time r org . The ergodic time is the character- 
istic time scale which determines how long it takes for an 
electron to reach boundaries of the system from arbitrary 
point in the grain, i.e. to explore the whole phase space of 
the system. The ergodic time is different in ballistic and 
diffusive systems. In ballistic systems, the ergodic time 
is determined by time between collisions with the system 
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boundaries, T org ~ L/v-p, on the other hand in the diffu- 
sive regime ergodic time is defined by the diffusion time 
between boundaries, T crg ~ L 2 /D, where D = v^lo/d is 
the diffusion constant. 

Another characteristic time of the system is the escape 
time t osc , which determines how long electrons spend in 
the grain before colliding with the contact. This time 
can be estimated as the ratio of the total boundary area 
and area of the contact, multiplied by the ergodic time. 
Indeed, the ergodic time determines time between col- 
lisions with boundaries, and the areas ratio determines 
the average number of collisions with closed boundaries 
before a collision with the contact takes place. We have 
t osc oc T elg (L/b) d ~ 1 . We quantify the escape time from 
the dot for a two dimensional system by considering the 
volume r of the phase space of the grain. The integral 
over the phase space of the dot can be calculated us- 
ing two different approaches. First, directly integrating 
over the direction of the momentum n and the volume of 
the grain, we obtain Y = 2nV g . On the other hand, ev- 
ery point in the phase space is uniquely determined by a 
trajectory, starting at the superconductor normal metal 
contact (SN contact). The integral over the phase space 
in this case can be represented as 



,7T/2 , 

r= / d6 cos6l{6,b)db = 2l esc S c 

and we obtain the escape time r osc = Z osc /i;f: 

1 v F S c 



Tesc — 



7rV g ' 



(1) 



(2) 



where V g is the area of two dimensional metal. Below we 
will usually use the escape rate 7 0S c = 1/t CS c- 

We emphasize, that the escape time t csc , Eq. (2), based 
on the trajectory calculation coincides with the escape 
time, defined in the random matrix theory as U/t csc = 
N5i/2tt, where N is the number of open channels in the 
contact, N = 26/Af, Si = 2Trh 2 /m e S g is the mean level 
spacing, Af is the Fermi wave length, m c is the electron 
effective mass and S g is the area of the dot. In this paper 
we assume S\ to be negligibly small and N to be large. 

The superconductivity order parameter A is constant 
inside the superconducting lead and vanishes fast at the 
contact with the metal grain. We consider the limit of 
strong superconductivity, A » 7 osc , and study the den- 
sity of states at energy e <~ 7 CSC much smaller than the 
energy gap A in the superconductor, £«A. 

The electron Hamiltonian is the sum of kinetic and 
potential terms: 



H = - £ F + V s (v) + V q (r), 
2m 



(3) 



where the potential term V s (r) + V q (r) is divided into 
semiclassical smooth V B part, which describes the effect 
of the grain boundaries, and the quantum part, V^, which 
represents impurity scattering: V q (r) = J2 i ViS(r — Hi). 



We write the full Hamiltonian in the Gorkov-Nambu 
space as 



Hqn = Ht z + A(r), 



(4) 



where the term A(r) = f x A#(r) takes into account su- 
perconductor pairing, f XjZ are the Pauli matrices and 
6>(r) = 1 inside the superconductor and vanishes other- 
wise. This approximation for the coordinate dependence 
of the order parameter is justified in the limit of strong 
superconductivity, A^> e. 

We introduce the Green's function as an inverse func- 
tion of the Hamiltonian: 

(e-H G ^G(s,v,v') = S(v-r')l (5) 

and perform the Wigner transformation 

G(e,R,p) = J e ipr G(e,R + r/2,R-r/2)dr. (6) 

Next, we change coordinates in the phase space from r 
and p to £ and x, where £ = H(r,p) — Ef and x = 
(r, n) with n = p/|p|. The reduced Green's function is 
obtained by integration over £: 



G(e,x) 



g{e,x) /+(e,x) 



- / G(e, 

TTP J 



(7) 

where v is the density of states in the grain per spin de- 
gree of freedom in the absence of superconductor. The 
semiclassical Green's function allows us to calculate elec- 
tron density of states in the system as an integral over 
the phase space of electrons in the grain: 



M{e) 



z/Rc J Tr{[f z ,g(e,x)]}^ 



(8) 



Here [•,•] stands for a commutator. An element of the 
phase space dx is dx = drdil n ; dil n denotes a measure 
on a sphere of a unit radius, representing electron mo- 
mentum direction n, and is the total area of a sphere 
of a unit radius in d— dimensional space, S% = 2ir and 
^ = 4tt. 

The reduced Green's function Q(e 7 x) satisfies the 
Eilenbcrgcr equation: 5,27 



ieT z +A(r),G(e,x) - CG(e,x) = I[Q{s,x% (9) 



where 



c = p d__ dv s (v) d 



m dr dr dp 



(10) 



is the Liouville operator and I[Q(e,x)] is the scattering 
term. The scattering term takes into account quantum 
part of the potential V q , introduced in Eq. (3). A general 
form for the scattering term in the Bohr approximation 



is: 



J[S(r,n)] - [E(r,n),0(r,n) 



(11a) 
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where S(r) stands for the integral 



S(e,r,n) 



C7 (£,r,n) <7+(e,r, n) 
cr_ (e,r,n) -er (e,r,n) 

Trimi Jg(r,n')\V q (n'-n)\ 



(lib) 



2 cZO n ' 



over direction of momentum n at the Fermi surface, 
Vi(n — n') represents the scattering amplitude off a single 
impurity, the impurity concentration is n\. 

A solution of the Eilcnberger equation, Eq. (9), is re- 
stricted by the constraints: 



g 2 (e,x) = l, tr£(e,x) = 0, 



(12) 



which represent the normalization conditions for Q(s, x). 

Inside the superconductor at distance £ away from the 
SN contact the Green function is given by 



</(£,£) = f x + a + e 



2A£ 



1 1 

-1 -1 



(13) 

where £ < for electrons moving in the superconduc- 
tor towards the grain and £ > for electrons moving in 
the superconductor away from the grain. The require- 
ment that the Green function inside the superconduc- 
tor relaxes to G{e,£) = f x , gives the following form of 
the Green function Q s (e,x) at the SN contact inside the 
superconductor: 18 



</ s (e,x) = t x + a 



1 



(14) 



where <; = +1 for electrons incoming to the grain, 
i.e. nn s > 0, and ? = —1 for outgoing electrons, 
nn s < 0. Here n s is a unit vector perpendicular to 
the superconductor-metal interface and directed to the 
grain. The second term of Eq. (14) has a specific form, 
so that it corresponds to a vanishing solution of the Eilcn- 
berger equation inside the superconductor. 18 Coefficient 
a is a free parameter. For a reflectionless SN contact the 
Green function is continuous and Eq. (14) represents the 
boundary conditions for electron Green functions inside 
the grain at the SN contact. 

We represent a point of the phase space by its coor- 
dinate r and the direction of the electron momentum n. 
Alternatively, we can parameterize the phase space by 
trajectories, which begin and end at the SN interface. 
Since only one trajectory goes through point r in di- 
rection n, the trajectory and the distance £ along this 
trajectory, measured from the starting point, uniquely 
represent point x = (r, n) of the phase space. Below we 
often use the trajectory parameterization of the phase. A 
point at the SN contact with nn s > has always C = 0, 
and a point at the contact with nn s < has £ = I, where 
I is the length of the corresponding trajectory. We omit 
the Fermi velocity and imply that the distance ( is mea- 
sured in time units. In this case, £ corresponds to time, 
which electron travels from the beginning of a trajectory 



TABLE I: The boundary conditions for the functions o±(e, £)■ 



C = o 



a+(e,C) 



a-(e,C) 



unrestrained 
1 



unrestrained 



to a given point, and I is the total time of motion along 
this trajectory. 

The normalization conditions, Eqs. (12), for the semi- 
classical Green's function C?(e, x) leaves only two indepen- 
dent functions, which parameterize the matrix Q(e,x). 
For certain purposes it is convenient to choose the Riccati 
parametrization in terms of functions a± = a±(e,x): 28 



1 



1 + a + a_ 



1 



- a+ci- 
2a- 



2a+ 
a + a_ — 1 



(15) 



The Eilenberger equations for functions o±(e,x) ac- 
quire the following form: 

(2ie — C)a + = o-a\ + 2a a + — cr + , (16a) 
(2ie + C) a_ = <j + a 2 _ + 2a a_ — er_. (16b) 

Here the self energy elements a± = cr±(e,x) and ero = 
CTo(e,x) arc defined by Eq. (lib). 

The boundary conditions, Eq. (14), acquire a simple 
form in the Riccati parameterization. Using the con- 
tinuity condition for the semiclassical Green's function 
G(e,x) at the SN interface, and comparing Eqs. (15) and 
(14) we find that a+(e, I) — 1 at the end of a trajectory 
of length I, [<; < in Eq. (14)] and a_(e,C = I) is not 
restrained and reflects a freedom of factor a in Eq. (14). 
Similarly, a_(e, 0) = 1 for incoming electrons, <r > 0, 
while o+(e, £) is a free parameter, see Table I. 

In the limit of strong superconductivity, A » e, an 
electron moves along some trajectory and when it is 
Andreev reflected as a hole, the hole moves in the op- 
posite direction along the same trajectory. Functions 
a±(e, r, n), which describe electron Green function in the 
Riccati parametrization, Eq. (15), have the meaning of 
the phase factor, which electron-hole pairs acquire as a re- 
sult of a single Andreev reflection at the SN contact while 
they travel from the contact to point {r, n} (note the 
boundary conditions in Table I). We refer to a±(s, r, n) 
as wave functions of an electron- hole pairs at point {r, n}. 
The off-diagonal elements of the Green function Q(e,x) 
may be represented as a series in a±(e, r, n) and corre- 
spond to the full amplitude, which takes into account 
multiple Andreev reflections, see Fig. 2. The diagonal 
element of the Green function may also be written in the 
form of a series, which gives us the following expression 
for the electron density of states in the grain: 

kn ^ 2 f dx 

x ( 1 + 2 E(- 1 ) fcRe { a +( e ' x ) a -( £ ' x )}) • 
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FIG. 2: In this figure we demonstrate the meaning of the 
Riccati parameterization, introduced by Eq. (15). Functions 
a± (e, x) represent a result of a single Andreev reflection on the 
SN interface. Particularly, o_(e, x) describes the process, in 
which an electron from point x travels to the SN contact and 
experiences an Andreev reflection, then a hole moves along 
exactly the same trajectory back to point x. In this figure 
electrons, represented by a solid line, move from right to left 
and holes, shown by a dashed lines, move from left to right. 
The elements /±(e,x) and <jr(£,x) may be represented in pow- 
ers of the a±(e, x) functions. 

Here 5i = vV g is the mean level spacing of electron states 
in the grain of the volume V g . The fcth term of this se- 
ries represents the effect of k Andreev reflections. We 
notice, that the first term of this expansion (unity) cor- 
respond to electron Green function in the grain without 
superconductivity at the contact. 

In conclusion of this Section, we make the following ob- 
servation. The scattering term vanishes after integration 
over the direction of momentum at each point of coordi- 
nate space. If we integrate the Eilenberger equation in 
the form of Eq. (9) over the phase space, we obtain 

/ nn s g(e,r,n)dS c ^- = ie [ [t z , Q{e, r, n)]dr^-, 

(18) 

where dS c is an element of the SN interface and n s is 
a unit vector, perpendicular to dS c . We notice that the 
diagonal part of Eq. (18) gives 

/ nn s g(s, r, n)dS c dfl n = 0. (19a) 
J s c 

The latter equation represents conservation of the elec- 
tric charge in the grain. The left hand side part of the 
off-diagonal part of Eq. (18) contains the following differ- 
ence f+(e, I) — f+(e, 0) of the off-diagonal elements of the 
Green function at the terminals of trajectory of length I. 
From the boundary conditions, defined by Eq. (14), we 
know that f+(s, Q and g{e,() ar e related at the trajec- 
tory terminals: 

/+(e,0) - 1-sM), 
/+(e,0 = l + g(e,l). 



These relations allow us to write the following equation 
for the off-diagonal part of Eq. (18): 

/ |nn s | g(s,r,n)dS c ^- = 2is [ f(s,r,n)dr^-. 

(19b) 

This equation makes a connection between the boundary 
values of the Green's function and the integral of the 
Green's function over the phase space in the metal grain. 

III. QUANTUM DISORDER 

In this section we study density of states in the dirty 
system, when the quantum part of the potential is pro- 
duced by small (point) impurities. The collision term can 
be written in the form of Eqs. (11). 

The mean free time is defined as the integral over dif- 
ferent directions of momentum: 

- = 2*vn, [ |I/(n-n')| 2 ^. (20) 

TO J 

Note, that the angular dependence of the Green's func- 
tion introduces scattering time, different from the mean 
free path To. Particularly, in the Usadel equation angu- 
lar dependence of the Green's function on the momentum 
direction leads to the appearance of the transport mean 
free time Tt r , see Eq. (22) below. 

In this section we assume that the scattering potential 
has no sharp angle dependence, and we use the mean free 
time as a characteristic value of the impurity scattering 
strength. 

A. Diffusive contact 

In the limit when the scattering is strong, Eilen- 
berger equation reduces to the Usadel equation 29 [we 
will discuss conditions for applicability of the Usadel 
equation below]. In this limit, the Eilenberger Green's 
function weakly depends on the direction of the mo- 
mentum n and we represent Green's function in the 
form: Q(e, r, n) = Oo(s, r) + ngi(e, r), where we assume 
Qo{e, r) 3> ngi(e, r). Taking into account the normaliza- 
tion condition, Eq. (12), we find 

gi(e,r) = -v F T tT g(e,r)-^-g(s,r), (21) 
where we introduced the transport mean free time r tr : 

— = 2um / (1 - nn')|U q (n' - n)| 2 — . (22) 

Ttr J "d 

Substituting gi(e, r) into the Eilenberger equation, 
Eq. (9), we obtain the Usadel equation for a system with 
time-reversal symmetry: 

£>V 2 e = 2is sine - 2A(r) cos 9, (23) 
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FIG. 3: The plot shows density of states in an Andreev 
billiard with strong disorder, when the mean free path lo is 
smaller than the size of the SN contact b [solid line]. The 
units of the horizontal axis e/7 represent energy in terms of 
7 = 7c see Eq. (27). For comparison, we plot the density 
of states obtained within a random matrix theory (see ref. 20 
and Sec. Ill B 1 of this paper) as a function of dimensionless 
energy 7 = 7 CSC , where 7 esc is the escape rate [dashed line]. 

where we parameterize the Green function by & = 9(e, r) 

a /cos9 sin6 \ , s 

G ={sme -cose ) (24) 

and D = VpTtr/d is the diffusion coefficient. 

We solve the Usadel equation for a d-dimensional metal 
grain of volume V g oc L d , connected to a superconductor 
by a small contact of the area S c , see Fig. 1. In the 
normal region we substitute A = into Eq. (23) and 
integrate the resulting equation over the volume V g of 
the grain: 

D ( n s V0(e,s) dS c = 2ie [ sin0(e,r)dr, (25) 

where the integral in the left hand side is taken over the 
area of the contact 5 C , and n s denotes the normal vector 
at the SN contact. Since 0(r) varies only close to the 
contact, we neglect variation of and replace the right 
hand side of Eq. (25) by 2eV g sineo, where ©o is an 
asymptotic value of 6(r) away from the contact. 

To calculate the left hand side of Eq. (25) we solve the 
following equation 

L>V 2 e = 0, (26) 

neglecting the term with energy e in the Usadel equation. 
This approximation is justified for energy smaller than 
the characteristic value of the left hand side of Eq. (26) , 
which we estimate as D/b 2 . 

We notice that Eq. (26) coincides with equation 
V 2 </?(r) — for an electrostatic potential ip(r) in the 
problem of calculation of the current 




through the contact. We can write the Ohms law / = 
(tpc — (po) / R c for the current /, where <p c and ipo are 
the values of the electric potential <p at the contact and 
inside the grain, respectively, and R c is the resistance 
of the contact. The contact resistance R c depends on 
the geometry of the contact. For a circular contact with 
diameter 6 we have aR c oc 1/6 for a three dimensional 
system and aR c oc logL/6 for a two dimensional system. 

Identifying the left hand side of Eq. (25) with the cur- 
rent I and the Usadel parameter 0(r) with the electric 
potential ip(r), we reduce the solution of the Usadel equa- 
tion to the electrostatic problem. Then, using the bound- 
ary condition for 6 = 7r/2 at the contact and introducing 
the value of the parameter e in the grain as © = tt/2 — #, 
we obtain the equation for 

tf = 2z-costf, lc = -R^. (27) 

Here we introduced the quantity j c , which has a mean- 
ing of the electron escape rate from the grain with a 
diffusive contact. Because a oc e 2 vD/h, we have 7 C oc 
(h/e 2 R c )(l/i/V g ) = iVtransCh, where TVtrans = h/e 2 R c is 
the number of transparent channels in the SN contact 
[compare j c to 7 CSC = N8\/2-n for a clean contact with 
N channels]. 

The density of states in the grain is given by the real 
part of the Eilenberger Green function 

M{e) = l-Resintf. (28) 

and vanishes, if d is purely imaginary. For small energy 
e < Eg , Eq. (27) has onlu imaginary solutions, and the 
spectrum begins at threshold energy Eg , where $ ac- 
quires a real part. From Eq. (27) we find the gap energy 

E^ ss 0.33l7c. (29a) 

Above the gap energy E^ the density of states exhibits 
a square root singularity: 

^ ^ il/iF 7 < 29b) 

and the numerical constant c\j ss 3.62. At higher en- 
ergy the density of states approaches its metal value u, 
according to the following expression: 

"<"-£( 1 + ss?)- <30) 

We plot the density of states, described by Eq. (28) 
through a solution of Eq. (27), in Fig. 3 by a solid line. 
For comparison, we present the random matrix result 20 
for the case of rcfiectionless channels by a solid line. 26 

This result is applicable for e <C f/6 2 . At larger ener- 
gies, e ~ -D/6 2 , Eq. (26) and, consequently, Eq. (27) for 
the parameter i9 arc no longer valid. 
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It is instructive to compare the result of Eqs. (29) with 
the solution of Usadel equation in the one dimensional 
case for wires of length L with both ends connected to a 
superconductor, see refs. [30,31]. The value of the energy 
gap and density of states above the gap are given by: 



1. 



(31) 



We observe that E g ~ D / L 2 . In this case the Laplasian 
term in Eq. (25) is comparable with the energy term and 
approximation of Eq. (26) is not applicable. The Green's 
function is obtained as a result of exact integration of the 
one-dimensional sine-Gordon equation. 30 ' 31 

Deriving Eq. (21) we have to assume that ng^ <C 1, on 
the other hand, using Eq. (21) and taking into account 
that V£(r) <~ 1/b, we obtain ngi <~ h/b. The result 
of Eqs. (29) and (30) is valid in the limit when l n <C b. 
Beyond this limit the Usadel equation is not applicable. 

For one dimensional systems, the density of states was 
studied beyond the Usadel equation in refs. [32], where 
an arbitrary relation between the system size L and the 
mean free path Iq was considered. The conclusion of 
refs. [32] is that there always exists a finite energy gap in 
the density of states. 



B. Ballistic contact 

In this subsection we study the density of states in 
the limit, when the mean free path is larger than the 
size of the contact Iq 3> b. In the intermediate regime, 
b <~ lo, the result depends on specific geometry of the 
contact, and the solution of the Eilcnberger equation may 
be studied numerically in the spirit of refs. [32]. 



1. Random Matrix Limit 

Here we consider the limit of short mean free time 
To *C t osc . The Green's function coincides with the Green 
function at the contact, see Eq. (14) only at small dis- 
tances C ^ To ° r I — C ^ T o along trajectories from the 
terminals of these trajectories. After several scatterings, 
at distances larger than the mean free path, the Green's 
function averages out and becomes independent from mo- 
mentum direction and coordinate. We refer to this coor- 
dinate independent part of the Green's function as a zero 
mode and define it by 



So 



9o(s) /o(e) 
/o(e) -ffo (e) 

1 f l-a 2 {e) 2a (e) 



l + a 2 (e) V 2a (e) -l + ag(e) 



• (32) 



Here the second equation introduces the Riccati param- 
eterization of the zero mode component in terms of a 
single function ao(e), see Eq. (15). 



The Green function G T (e, () near the boundaries (C = 
0,0 is: 



G(e 0) - G T (e I) - ( ff(£ ' 0) /+(e ' 0) 

= 1 ( 1 - a (e) 2a (e) 
l + a (e)V 2 -l + ao(e) 



(33) 



Here we took into account the boundary conditions 
a_(e) = 1 for incoming electrons and a + (s) = 1 for out- 
going electrons, see Table I. The unrestrained values of 
a + (e) and a_(e) at the SN contact coincide with the zero 
mode value a (e) inside the grain. 

At sufficiently small values of the mean free time t , 
we neglect the contribution to the average value of the 
Green's function from the parts of trajectories close to the 
SN contact, where G(e, () is given by Eq. (33), since this 
contribution is as small as 7 CS cTo, where 7 0SC is defined 
by Eq. (2) in terms of the contact area and the grain 
volume. Then, the integral equation (19b) reduces to 



7esc3(e,0) = ief (e). 



(34) 



Substituting /o(e) from Eq. (32) and g(e, 0) from 
Eq. (33) into Eq. (34), we rewrite Eq. (34) in the form: 



1 - a (e) 



ie- 



2a (e) 



(35) 



l + a (e) "~l + al(s)' 
We substitute a (e) = cxp(— 2i^) and reduce Eq. (35) to 



tan ip cos 2tp = . 

7esc 



(36) 



The density of states 



Af(e) = |-Re[(?o(e)] = |"Imtan2^ (37) 

is defined by Eq. (32) in terms of the solution a (e) and 
is finite only if ip has an imaginary part. The function in 
the left hand side of Eq. (36) has a maximum at certain 

value Vm = arccos(\/l + v / 5/2) ~ 0.45 and we define 



tan ipm cos 2V>m7c 



'5\/5- 11 



7csc w 0.30 7cS c. (38a) 



At small energy e < E™* , Eq. (36) has real solutions and 
the density of states vanishes. Above the gap we expand 
the left hand side of Eq. (36) in tp — ip m as v(-p — ip m ) 2 = 
i?™* — e, with v ~ 2.17, and obtain that the density of 
states has a square root singularity above the gap 



i 



(38b) 



with c rmt = 2 J 2 + 4/V5 w 3.89. 
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At energies larger than the escape rate 7 eS c, the den- 
sity of states approaches the normal metal value 2/(5!, 
according to 



Af(e) = - 

Ol 



1 + 



7 2 

I CSC 

2e 2 



(39) 



The plot of the density of states is presented in Fig. 3 by 
the dashed line. We emphasize that the above result for 
the density of states, Eqs. (36)-(39) coincides with the re- 
sult of rcf. [20] obtained within random matrix theory. 21 



2. Crossover regime 

Here we consider arbitrary mean free time r ~ r osc 
due to scattering off isotropic impurities. In this case we 
have the following form of the self energy: 



(40) 



i.e. the self energy is proportional to the angle average of 
the Green's function at point r inside the grain. For small 
dots, when the Thouless energy F/ph is the largest energy 
scale of the grain, we will make a conjecture, that the 
averaging over direction of electron momentum is equiv- 
alent to averaging over the full phase space of the system: 




esc 



FIG. 4: The density of states for different values of mean 
free time to is presented. In dirty metal grain (To7 esc <C 1) 
the density of states is described by a random matrix theory 
[solid line]. As the strength of disorder decreases, the gap size 
decreases and the rise of the density of states smoothes. The 
cases of moderate disorder, To7esc = 0.7, and weak disorder, 
To7esc = 10, are shown by a dotted and dash-dotted lines 
respectively. [The density of states for ro7 eS c = 0.7 is shown 
only near the gap.] Only in the semiclassical limit the density 
of states has no gap, although it is exponentially suppressed 
at small energy. 



<S(e,r,n)) 



1 



Q(e,r, n)dr 



/ 



= G(e,r,n) 



(41) 



Indeed, different directions of the momentum represent 
contributions from different points of the grain, since we 
assume that the Green's function does not change fast 
at short parts of trajectories, corresponding to electron 
motion during the ergodic time r org = l/£rh- 

The Green's function <j(e,x) satisfies the Eilenberger 
equation in the form: 



CG(e, r, n) = [H cS (e), G(e, r, n)], (42a) 



where 



^=( g fI X) =i ^ + k {C{£ ^ ]) - (42b) 

Assuming that H c s{e) is fixed, Eq. (42a) is linear in 
Q(e, r, n). We choose the parametrization of the grain's 
phase space in terms of the classical trajectories, which go 
through space point r and in the direction of n. Then, 
the solution of Eq. (42a) can be found for each trajec- 
tory as a function of the coordinate along the trajectory 
(, since the Liouville operator has the form C = d/dQ. 
In this trajectory parametrization of the phase space, the 
smooth part of the potential V s is taken into account as 
the refraction of trajectories. 



The solution to Eq. (42a) can be written in terms of 
three matrices: 



-f-Hett] Q± 

L>0 



1 



2A, 



—Fq Go ± D 
Go T -Do Fq 



.(43) 



Matrices g z and g± are related to the Pauli matrices f z 
and t± = (f x ± ir y )/2 through a rotation in the Gorkov- 
Nambu space, which diagonalizes the Hamiltonian H c s 1 
see Eq. (42b). If H e s is diagonal, matrices Qi coincide 
with matrices fj . A general solution of Eq. (42a) is 

Q(s, C) - a z g z + a +Q+ e 2D ^- 1 ^ + a ^^- 2D ^- l ' 2 \ 

(44) 

Fq and coordinate £ is a length 
along a trajectory of length I. Coefficients a z and a± 
are uniquely determined by the boundary conditions, 
Eq. (14). We have 



where Dq 



r 2 
u o 



OL z = 



a± = 



F cosh D l + D sinh D l 
D cosh D n l + Fo sinh D l ' 

Go 

D cosh D l + Fo sinh D l 



(45a) 
(45b) 



The average of the Green's function is defined as the 
integral over the phase space, which in the trajectory 
parameterization of the phase space is 

/•oo A 

0(e,r,n)) = 7esc / diP(l) &{e,x)dC, (46) 
Jo Jo 
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where we first integrate Q(s,x) along a trajectory of 
length I, and then average the result over trajectories 
of different length I with the weight 



P(l) = 7escexp(-7e S cO- 



(47) 



Here 7 esc is escape rate, or inverse average escape time 
t osc , defined by Eq. (2). 

Performing the integrations, we obtain: 



1 



?7i sinh 9 — i]2 i]i + i]2 sinh 9 
cosh 9 \ m + m sinh 9 r] 2 - rji sinh 9 



and coefficients 771.2 are 



(48) 



m 



m = 



+ 



1 + 2E 



(-l) fe 7, 



k„,2 

CSC 



—J (2F kcosh9 + 7 CSC ) : 



tanh 2fc -, 



k 

7csc , - tanh- 
F coshfl 2 

9 00 
■^7esc \ 



(-l) fe 7c 



F sinh 8f^ 2F k cosh (2 + 

/C; 

/C— 1 



■ tanh 



2fc 



Here we introduced a shorthand notation 9, defined 
through Go = F sinh# and D = F cosh 6*. 

In these calculations, Eqs. (42b) and (48) represent a 
self consistent system of equation for Fq and Go (or 6*). 
The solution of this self consistency equation gives the 
electron density of states in the grain: 



Af(e) = —Re \ i]i tanh( 



1 



cosh0 



(50) 



The density of states, in the form of Eq. (50), sup- 
plied with Eqs. (42) and (48), are obtained here for the 
first time. These equations can be solved numerically for 
an arbitrary value of To7 OS c- In Fig. 4 we plot M(e) for 
weak disorder, r 7csc = 10, and for moderately strong 
disorder, To7 eS c = 0.7, by dash-dotted and dotted lines, 
respectively. For weak disorder, the density of states is 
similar to the semiclassical result [dashed line] , although 
it exhibits the gap E g ss 0.097 esc . The gap is not notice- 
able in Fig. 4, since the overall number of states near the 
gap is already exponentially small. As the strength of 



disorder increases, r n 1 > 



7e 



the density of states ac- 



quire the form, similar to the random matrix result [solid 
line], see the dotted line in Fig. 4 for the To7 OS c = 0.7 case 
shown only for e < 0.57 OSC . We also present the energy 
gap E g as a function of the mean free time To. As To in- 
creases, the gap size monotonically decreases, see Fig. 5. 

Now we show, that both the semiclassical result 18 ' 19 
and the random matrix result 20 may be obtained from 
Eqs. (42), (48) and (50). 

In the absence of impurity scattering, To — > 00, we 
obtain the ballistic semiclassical result. Indeed, we have 
Go = -Do = is, Fo — * and 9 — > 00. The density of states 
is determined by 771: 



(-l) fe 7, 



2 

esc 



2 — 

^^T,^ (27fe + 7csc) 2 ' 

k= — 00 



(51) 



0.3 



w 



0.2 



0.1 



0.0 



T Y 
0'esc 



FIG. 5: The gap energy E s as a function of the mean free 
time to is shown by a solid line. 



This result is equivalent to the semiclassical expression 
for the density of states in the Andreev grain: 



A/Ye) = 2 7r2 7csc cosh7r7e SC /2g 
1 ' 8 X 4e 2 sinh 2 7T7 csc /2 ; r' 



(52) 



which was derived in ref. [19]. At large energy, the density 
of states in the clean semiclassical system approaches its 
value in the normal state according to 



(53) 



In the opposite limit, er <g; 1 and 7 osc To <C 1, we 
neglect series in Eq. (42b) for 7,72, since the series are of 
the second order in 708^0- To the zeroth order in etq and 
7oscTo the self-consistency equation, Eq. (42b), is satisfied 
by arbitrary 9. Keeping the first order terms in £To and 
7 esc To, we obtain the equation for 9: 



is — 7csc c°sh 9 tanh • 



and the expression for the density of states 



A/» = y- Re {tanh . 
01 



(54) 



(55) 



Equations (54) and (55) are equivalent to Eqs. (36) and 
(37) after the substitution 9 = 2iip. We conclude, that 
the density of states to the lowest order in et and 
7cscTo corresponds to the random matrix theory result, 
see Eqs. (38) and (39) and ref. [20]. 

The most drastic difference between the density of 
states corresponding to different values of the mean free 
time To appears at small energy, see Figs. 4 and 5. This 
difference was discussed in refs. [18,19,20,22,23,24]. Here 
we emphasize that the difference exists even at high en- 
ergy. Comparing Eqs. (39) and (53), we notice that al- 
ready the lowest order in 7 0S c/e term differs by a numer- 
ical factor of 7r 2 /12. In general, a high energy asymptote 
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£ S> 7csc of the density of states can be written in the 
form: 



2 



7 2 

\ / CSC 



where (p(er n ) is 



f{y) = R-c 



1 + 2i - 2y - ttV/12 

(i + iy) 2 



(56a) 



(56b) 



We observe, that transition from the dirty limit, etq <C 1, 
to the clean limit, er S> 1, occurs at energy of the order 
of the scattering rate, e <~ l/r . 

We discuss the difference between the high energy lim- 
its of the semiclassical and the random matrix theory 21 
results in more details. Mathematically, the different co- 
efficient in front of £ 2 /7csc m Eqs. (39) and (53) may 
be attributed to J2(~K ) l k2 

over k, when the random 
matrix theory result is obtained only from the first term 
of the sum, while the semiclassical result corresponds to 
X;^Li(-l) fc /fc 2 = -7r 2 /12. We argue that this mathe- 
matical difference between Eqs. (39) and (53) has a phys- 
ical meaning. 

In the semiclassical limit, the energy levels in the grain 
are usually interpreted as standing waves along classi- 
cal trajectories. 18 ' 19 Alternatively the contribution to the 
density of states may be explained as the effect of inter- 
ference of repetitive motion of electron-hole pairs along 
a classical trajectory, see Eq. (51). A scattering off im- 
purities switches electrons and holes between different 
trajectories and destroys the interference effect of such 
repetitive motion. Due to the scattering, an electron- 
hole pair travels along any trajectory only once, and the 
sum over k in Eq. (51) - the number of repetitions - is 
limited to k — ±1. 



IV. QUANTUM CHAOS 

In this Section we consider a clean metal grain with 
chaotic smooth potential and consider the effect of quan- 
tum diffraction on electron density of states. 



A. Qualitative discussion 

In quantum mechanics electron propagation is de- 
scribed by the evolution of a wave packet. The electron 
Green function, see Eq. (7), at some point x of the phase 
space is actually determined by averaging over a set of 
wave functions, which constitute a wave packet with cen- 
ter of mass at point x. In this sense, the effect of quan- 
tum diffraction is similar to the effect of disorder. We 
emphasize, however, that the quantum diffraction mixes 
wave functions of electron- hole pairs at close points of the 
phase space within a minimal wave packet, see below. On 
the contrary, scattering off point impurities in disordered 
system mixes electron-hole wave functions propagating 




b) 



FIG. 6: Figure a) represents a\{e, x)a 2 _ (e, x) for the case, 
when the interference exists. The low figure corresponds to 
the in dependent propagation of two electron-hole pairs. In- 
deed, if one pair is Andreev reflected at point D, the other is 
normal reflected at point B and the paths, along which these 
two pairs travel from points B' and D' , are different by a 
length of order of l CBC . The paths AA' , BB', CC and DD' 
are straitened and the electron-hole pairs have multiple re- 
flections on the dot normal boundaries while they are moving 
along these paths. 



in completely different directions. Thus, the scattering 
off point impurities results in the mixing of distant points 
of the phase space. 

First we estimate the size of a typical wave packet in 
the phase space and then discuss the effect of chaotic 
dynamics on the electron Green function. 

We consider a wave packet, which consists of two 
electron-hole pairs, travelling from points C and A' to 
points B 1 and D' in Fig. 6. The width of the wave packet 
is Sp and pairs' momenta are misaligned within angle Scj), 
so that the momentum distribution of the wave packet 
has width p± = pf5<P, where pf = 2ttH/\f is the Fermi 
momentum and Ap is the Fermi wavelength. According 
to the uncertainty principle, we have p±Sp <~ 2nh. Tak- 
ing into account the relation Sp s=a LScj), see Fig. 6, we 
obtain 



= y/X F /L, 



(57) 



for the momentum deviation angle and for the coordinate 
width of a minimal wave packet, respectively. 

The Eilenberger equation, Eq. (9), for a clean system 
does not describe quantum diffraction. The way to mimic 
quantum diffraction is to introduce small angle scattering 
potential, see refs. [7,9]. This scattering introduces mix- 
ing, or "interaction", between electron-hole wave func- 
tions, which approach each other within a small distance 
in the phase space, described by Eqs. (57). The small 
angle scattering suppresses the oscillating components of 
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the Green functions and destroys the semiclassical de- 
scription of electron-hole motion in the system. 

We develop the application of the method of small an- 
gle potential to calculations of density of states in An- 
dreev billiard in Appendix A, using the technique de- 
veloped in ref. [7]. Below in the text we utilize some- 
what more intuitive approach. We represent the Green 
function in the Riccati parameterization, see Eq. (15) 
and rewrite it in terms of well defined series in a±(e, x). 
To include the effect of quantum diffraction, we identify 
powers of a±(e, x) as weighted average over the minimal 
wave packet of the electron-hole pairs, see Fig. 6: 

k 

a±(e,x) -> (a±(e,x)) = Y[ / ^(x- x l )a ± (e,x l )dx 4 . 

(58) 

Function P(x — x') represents the coordinate and mo- 
mentum distribution within a minimal wave packet and 
vanishes fast on coordinate scale of p m in or for momentum 
deviation on angle larger than <5</> m in, see Eqs. (57). A 
detailed form of the function P(x — x') is not important, 
since, as we show below, the density of states depends 
logarithmically on the width of the function P(x — x'). 

In chaotic system the distance between two trajectories 
grows exponentially in time t: p(t) = p m i n exp(Ai) and 
4>(t) = 4> m i n exp(At), where A > is the Lyapunov expo- 
nent of the classical chaotic motion. Consequently, even 
though the size of the wave packet is initially small, time 
evolution leads to the spreading of the wave packet over 
the phase space. Generally, the spreading of the wave 
packet is characterized by the Ehrenfest time ts- The 
Ehrcnfest time shows how long it takes for a wave packet 
to spread over the whole phase space of the system and 
is equal to 

t E = \\n^ = hn—. (59) 

A Ap A Pmin 

We argue, however, that the relevant time scale for the 
problem of proximity effect in Andreev billiards depends 
on the size of the contact b and is given by 

T E = -ln = ^--T ln T' ( 60 ) 

A Pmin 2 A b 

The second term in the right hand side of Eq. (60) 
may be omitted, since it contains a logarithm of macro- 
scopic parameters b and L. Nevertheless we explain the 
appearance of the scale b as follows. If electron prop- 
agates from points B' or D' along trajectories B'B or 
D'D, see Fig. 6a, for time smaller than te, the spread- 
ing of the initial wave packet is smaller than the size 6 
of the SN contact and the two trajectories ends nearly 
simultaneously at points B and D. We conclude that 
the contribution of electron-hole wave functions to the 
electron Green function originates from coherent motion 
of electron-hole pairs along these two trajectories. This 
coherent contribution may be viewed as a repetitive mo- 
tion along the same trajectory, see the text at the end of 
Sec. IIIB2. 



On the other hand, if electron propagates along B'B 
and D' D for time longer than te, the wave packet spreads 
on the scale larger than the size b of the SN contact. As 
shown in Fig. 6b, although the electron-hole pair travel- 
ling along D'D experiences Andreev reflection at point 
D, the other pair travels along B'B and misses the SN 
contact. It continues its motion further from point B and 
explores the phase space, different from the one reached 
by the D'D trajectory. Thus, even though the wave 
packet originates from the contribution of electron-hole 
wave functions at close points in the phase space, the cor- 
responding wave functions are quite different, since they 
represent motion in distant parts of the phase space. The 
averaging over the wave packet leads to appearance of a 
component of the Green function, produced from differ- 
ent parts of the system phase space, known as a "zero 
mode" . 

We expect that if the Ehrcnfest time r E is much longer 
than the time of electron propagation along a typical tra- 
jectory, then the electron Green function may be deter- 
mined by semiclassical methods, see refs. [18,19], while 
in the opposite limit of short te, the quantum diffrac- 
tion destroys the repetitive modes and leads to strong 
coupling between electron-hole wave functions from dis- 
tant parts of the phase space. In this sense, motion of 
electron-hole pairs in systems with short te «C r esc actu- 
ally coincides with motion in disordered systems, and the 
random matrix result 20 describes the density of states, 
see Sec. IIIB1. 



B. Quantitative analysis 

In this Subsection we perform an analytical analysis of 
the density of states in clean chaotic Andreev billiards. 
First we consider the high energy limit e » 7 0SC , when the 
perturbation theory is applicable. In the perturbation 
theory we identify the relevant features of the Eilenberger 
equation for a description of motion of electron-hole pairs 
in the regime of quantum chaos. Next, we construct a 
self-consistency equation, which determines the density 
of states at arbitrary energy e <~ 7 osc . 



1. High energy limit 

As we showed in section IIIB2, the density of states 
of disordered systems approaches the normal density of 
states as 7e SC / e2 at large energy e » 7 0SC , but the coeffi- 
cient of this asymptote depends on the strength of disor- 
der, characterized by the mean free time To, see Eq. (56b). 
A similar asymptotic behavior exists in quantum chaotic 
systems and now the coefficient depends on the value of 
the Ehrcnfest time te- To prove this statement, we cal- 
culate the density of states at high energy (e » 7 0S c)- 

We use the Riccati parameterization of the Green func- 
tion and apply Eq. (58) for calculation of powers a±(e, x) 
at coincident points, which appear in the expression for 
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the density of states, see Eq. (17). A more formal deriva- 
tion with Perron-Frobenius differential operator is pre- 
sented in Appendix A. 

The density of states is determined by the sum of 
(a+(£, x))(a?l (e, x)) over all positive integer k. Here (. . . ) 
stands for the averaging over the phase space of electrons 
in the grain, while (...) denotes averaging within a min- 
imal wave packet according to Eq. (58) . We observe that 
a\{e, x) and a k _(e, x) are mutually independent, since the 
corresponding electron-hole pairs move in completely dif- 
ferent parts of the phase space, approaching point x from 
two opposite directions (notice the different signs before 
the Liouville operators in Eqs. (16)). Due to indepen- 
dence of a+(e, x) and a^_(e, x), we write 

(a* (e,x))(a* (e,x)) = (a* (e,x))(a* (e,x)), (61) 

so that the density of states in the grain is determined 
by the product of (a+(e, x)) and (a k L(e, x)). 

We start with cal culation of the averages of electron- 
hole wave functions a±(e, x). We integrate both sides of 
Eq. (16) over the phase space of electrons in the grain 
Vgfld- As we already mentioned, the scattering term 
vanishes after integration, see Eq. (18). The integral 
of £ cm a±(e,x) reduces to the values of a±(e,x) at the 
SN contact. At e ^> 7 OS c, only the contribution from 
the trajectory terminal with restrained value of a±(e, £) 
is important, since the unrestrained value of a±(s,Q at 
the SN contact oscillates fast and thus averages out. We 
obtain 

Mi^0«^- (62) 

Next, we calculate the contribution to the density of 
states from the k — 2 term in Eq. (17), (a\(e, x)), to the 
second order in 7 CS c/£- We show that this term behaves 
differently in semiclassical and quantum mechanics. Par- 
ticularly, the quantum diffraction destroys the semiclassi- 
cal contribution to the density of states, originating from 
the repetitive motion. 

A similar problem was solved in ref . [7] , where the weak 
localization correction to the conductance of quantum 
chaotic systems was studied. It was shown that the weak 
localization correction originates as a result of quantum 
diffraction which couples Cooperon and diffuson modes. 
Following the technique, developed in refs. [7,9], we in- 
troduce new functions A^(e,xi,x 2 ) with arbitrary sepa- 
rated points xi = (ri,ni) and x 2 = (r 2 ,n 2 ): 

Af(e,x 1 ,x 2 ) = o±(e,xi)a±(e,x 2 ). (63) 

We arc interested in the limit of small relative angle 4> 
between vectors ni and n 2 within a minimal wave packet, 
Eq. (57). We change variables from Xi and x 2 to the co- 
ordinates of the center of mass R = (ri + r 2 )/2, N w 
(ni + n 2 )/2 (for \cf)\ <C 1) and to the relative coordinates 
p = Yi — r 2 and <f). Under the condition \(j)\ <C 1, the 
vector of relative distance p is perpendicular to the di- 
rection of N and we refer to it through its projection p 



on the axis perpendicular to N. At \(j>\ <C 1 the equation 
for A 2 t (e,xi,x 2 ) in new variables is 

(4ieT£ cm -jC(p,<P))A±(e,R,-N,p,<f>)=0, (64) 

where 

r d 1 0U S d 

C cm = v F N— — — (65a) 

aR pf aR oN 

is the part of the Liouville operator describing the motion 
of the center of mass, and 

d 1 d 2 U d 
£(^)^---^- (65b) 

is the Liouville operator of relative motion of two 
electron-hole pairs. 

We average Af(s, R, N, p, 4>) over the position of cen- 
ter of mass N and R in the phase space. Integration of 
the second term with the Liouville operator £ cm by parts 
gives values of A^ at the boundaries. Particularly, due to 
the boundary conditions, see Table I, for outgoing elec- 
trons a±(e, ( = I) — 1 and we have A^e, R, N, p, <fi) = 1 
at the SN contact, provided Nn s < (n s is a normal vec- 
tor to the SN inerface). In the opposite limit Nn s > 0, 
the value of ^4^(e, R, N, p, <f>) is unrestrained. Similar 
analysis provides values of A^ie, R, N, p, <j>) at the SN 

contact. We introduce notation s 2 = (Ao) SIL /Af for the 
ratio of the average unrestricted value Af at the SN con- 
tact, (A^) sn , to the average of Af over the whole phase 
space. We obtain 

Lie+^TCM)) A$(e,p,<t>) = — . (66) 

At high energy s 2 w 1. 

The Liouville operator £(p, <j>) in Eq. (66) can be 
rewritten in new variables z and \-> introduced as 

1, <f> 2 L 2 + p 2 <t>L 
z = -ln — , x = arctan— . (67) 

As was shown in ref. [7], the variable \ is irrelevant for 
quantity Af(s, p, cf>), described by a solution of Eq. (66). 
Indeed, we represent evolution of the wave packet in time 
t, which is a conjugate variable to energy variable e in 
the left hand side of Eq. (66), see Fig. 7. The character- 
istic scale of the wave packet, represented by the variable 
exp(z) grows in time as exp(Ai) with the Lyapunov ex- 
ponent A. A particular configuration of the wave packet 
is described by variable \, which significantly fluctuates 
as the wave packet explores the phase space and gets 
averaged out after integration of A^ over a position of 
the center of mass. [Actually this averaging occurs at 
crgodic time scale, which is much shorter than any other 
time scale in the problem.] The fluctuations of the Lya- 
punov exponent A may also be taken into account. The 
resulting form of Eq. (66) is 7 : 

( 4ie + 7 csc - A— - Af(e,z) = 7 CSC , (68) 



13 



H 




^ >^ exp(z) 














J 'l4> 

















l/A 


2/A 


3/A * 



FIG. 7: Evolution of the wave packet in time. The ellipse 
represents the wave packet in the phase space. The variable 
2 determines the size of the ellipse (its diagonal), while the 
angle variable \ determines the shape of the ellipse. The 
dashed line in t-p plane represents the Lyapunov divergence 
of two trajectories p(t) = p m i n exp(Ai) with the Lyapunov 
constant A. 



where A is the average value of the Lyapunov exponent 
over the electron phase space in the grain and A 2 rep- 
resents the deviation of Lyapunov exponent for different 
points of the phase space, see ref. [7] for more details. 

As the size of the wave packet reaches the size of the 
SN contact b, and z ~ Inb/L, the a±(e, x) functions in 
become independent. 



A 2 (e, xi,x 2 ) become independent. Indeed, these func- 
tions are solutions of the Eilenberger equation along dif- 
ferent trajectories of significantly varying lengths. We 
may write the following boundary condition: 



Af (e, z = In b/L) = a± (e, x) 



(69) 



where a±(e,x) is the average of a±(e,x) over the phase 
space and is given by Eq. (62). [We should distinguish 
the average o±(e, x) and the zero mode components of 
the Green function. More precisely, the right hand side 
of Eq. (69) contains the zero mode value of a±(e, x), as 
we will discuss in the next subsection. At high energy 
£ S> 7esc both quantities are equal.] 

To solve Eqs. (66) and (69), we represent Af(e,p,<j>) 
as a sum of two functions: 



At(e,p,<f>)=A± b (e)+A± q (e,z), (70) 



where 



At h {e) = 



lose 



4«e + 7 CSC Aie 

is a solution of Eq. (68) without any boundary condi- 
tions imposed and A^ q (e, z) is a solution of homogeneous 
equation Eq. (68) with certain boundary conditions at 
z = Inb/L, so that A^ie, z) in the form of Eq. (70) sat- 
isfies Eq. (69). 



A general solution of the homogeneous Eq. (68) was 
found in ref. [7] in the limit 7 CSC , s <C A: 



w(lo, z) = exp 



ILU 

T 



X 2 uJ 2 



(72) 



with LO = 4e — 27 CS c- Taking into account the boundary 
conditions, we find 



AUs,z) 



(a±(e,x) 2 -^ 6 (e)) 



x w{Ae - i~f CSCl z - Inb/L). (73) 

This equation describes the Lyapunov divergence of 
classical trajectories, once a finite displacement between 
two trajectories was created by quantum diffraction. As 
was discussed in ref. [7], Eq. (72) neglects the effect of 
quantum diffraction and a more accurate solution for 
the function w(ui, z) is needed. 7 Indeed, according to 
Eq. (17), the density of states contains term A^(e,p,<f>), 
taken at coincident points, i.e. p, <j> — > 0. This limit corre- 
sponds to z — > — oo, when w(e, z — ► — oo) = 0. The latter 
equation means that the effect of quantum diffraction, 
which couples Green functions at different semiclassical 
trajectories is not described by the Eilenberger equation 
without scattering term (l/r q = 0). The scattering term 
modifies the equation for function w(e, z), see Appendix 
A for details. We argue that instead of working in the 
limit z — > — oo, within a logarithmic accuracy we can av- 
erage quantity Af(e, p, <f>) over p and <p within a minimal 
wave packet, Eq. (57). We write 



r 2 (e) = {w{4e-ij esc ,z-ln-jj. (74) 
Performing averaging over the minimal wave packet 
{w(u,z)) - J 'w(w,z) e-^ 2L2 +P 2 y 2 P 2 ^d(j)dp 



with z = In ((f? L 2 + p 2 )/b 2 we obtain 
r 2 (e) = exp -7o SC te - 4z£te 



8A 2 e 2 



A 2 



-TE 



(75) 



and te = (l/A) ln6/p m i n is the Ehrenfest time, defined 
by Eq. (60). Combining Eqs. (62), (71) and (73) we find 

(^( £ ^)) = ^{i-r 2 ( £ )} (76) 

to the lowest order in 7 csc /£. 

Equation (76) has a simple interpretation. Function 
T 2 (e) is the probability amplitude, that two electron-hole 
pairs, initially situated within the same minimal wave 
packet at point x of the phase space, escape the grain 
independently. In other words, r 2 (e) is the probability 
amplitude that a minimal wave packet, constructed at 
point x expand to size b of the contact before it leaves 
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the grain. Consequently, 1 — I^fe) is the probability am- 
plitude for the wave packet to leave the grain before it 
decays due to quantum diffraction. 

We generalize the above calculations of (A^_ 2 ) for ar- 
bitrary integer k > 1. We define 33 

k 

A±(e,{^}) = l[a ± (e,x l ) (77) 



and write the equation for (e, {xj}), averaged over the 
center of mass of variables {xi}: 

( 2 kie+ — T£(p,<t>))A± (s,p, <(>) = — , (78) 

\ 7"csc / 7"csc 

which is a generalization of Eq. (66) to arbitrary k. Here 
the variables p = {pi} and 4> — {4>i\ constitute a set 
of 2{k — 1) variables. We argue that 2k — 3 of these 
variables determine the relative position of k electron- 
hole pairs within the corresponding wave packet and thus 
are irrelevant, while one variable determines the scale of 
the wave packet in the phase space z = h\(L 2 ^4>1 + 



pf)/2, see Fig. 7. When the wave packet reaches the 
size of the SN contact, the wave functions of electron- 
hole pairs a±(e 7 x^ become independent and we have the 
boundary condition z = lnfe/L: 



At(e,z = lnb/L) = a ± (e, X ) cx -^L. (79) 

(IS) 

Following further along the steps of calculations of 
(Af(e, z — lnfe/L)), we obtain to the lowest order in 

7csc/e 

( ^ ±(e ' z)> = ii {1 " rfe(e)} (80) 

with 

/ 2fc 2 A 2 £ 2 \ 
T k (e) = exp ( -7c SC t-e ~ 2iker E — — r E J . (81) 

Using Eq. (17) we obtain the density of states at large 
energy, e » 7 0SC in the form 



Af(e) 



2 




OG 

k=2 



(-1)* 

k 2 



Re < 1 — exp 



2ike + 7c 




(82) 



We found that the averages of functions (a+(e, r, n)) 
are expressed in terms of Tk(s). Particularly, function 
r 2 (e) also appears in the expressions for i) the weak lo- 
calization correction to the conductivity 7 , ii) level statis- 
tics correlation function 8 , Hi) noise of a current through 
a chaotic cavity. 9 

Functions Tk(e) oscillate with a period proportional 
to Tg 1 . Similar oscillations were found in the weak lo- 
calization correction to the conductivity and in the en- 
ergy level statistics. 8 The origin of these oscillations has 
the following explanation. The effect of quantum diffrac- 
tion can be neglected for motion at time shorter than 
the Ehrenfest time. Nevertheless the diffraction effects 
appear nearly instantaneously, with characteristic turn- 
on time equal to A -1 , as time of motion approaches the 
Ehrenfest time. Thus we can separate two types of con- 
tributions to the Green function. One contribution rep- 
resents semiclassical motion of electron-hole pairs and 
is responsible for oscillations in the density of states in 
the phase space. The other contribution is a constant 
in phase space component of the Green function, which 
corresponds to the zero mode. 

To illustrate the oscillations of the density of states, we 
plot the tail of the density of states Af(e) at large energies 
£ S> 7csc for te = 1.5t osc , shown by a solid line in Fig. 8. 
For comparison, we also plot the density of states for the 



random matrix limit (dashed line) and semiclassical limit 
(dash-dotted line). 

We emphasize that it is the sharp turn-on of the ef- 
fects of quantum diffraction at the Ehrenfest time leads 
to oscillations of the density of states. Fluctuations of 
the Lyapunov exponent suppress these oscillations. In- 
deed, the Lyapunov exponent varies for different points 
of the phase space and its fluctuations affect functions 
Tk(e) through the term A 2 £ 2 te/A 2 . Thus for small ener- 
gies, e <C a/ A 2 /A2Te, we can disregard fluctuations of the 
Ehrenfest time. As energy increases, the oscillations of 
the density of states become suppressed, see Eq. (82). We 
also notice that in the quantum disorder, the process of 
scattering off impurities is Poissonian, and fluctuations 
of the free path of the order of the mean free path it- 
self. Consequently, oscillations of the density of states in 
disordered systems are absent, see Eqs. (56). 



2. Arbitrary energy 

Although the high energy expansion cannot be applied 
to the calculation of the density of states at energy com- 
parable with the escape rate 7 esc , exactly the same mech- 
anism of quantum diffraction is responsible for the ap- 
pearance of the gap in the density of states. 
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FIG. 8: The tail of the density of states at high energy, 
£ S> 7csc- The density of states, described by Eq. (82) for te = 
1.5T eS c, is shown by a solid line. The semiclassical density of 
states, te — > oo is shown by dash-dotted line, and the random 
matrix theory result, te — > 0, is presented by a dashed line. 
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FIG. 9: The dynamics of a minimal wave packet, initially 
situated at point P or P' is describe by the characteristics 
(dashed lines) of Eq. (83) are presented in £ — z. The bound- 
ary conditions are defined at the terminal of the character- 
istics. If characteristic intersect z — \nb/L at £ < I, we say 
that the electron-hole wave functions at the intersection are 
described by a zero mode solution, A~^(e, £, In b/L) — a k (e), 
while if the characteristics ends at the SN contact, the bound- 
ary condition is A£ (e, I, z) — 1. 



Equation for Af (p, cj>) becomes significantly simplified 
after averaging over the phase space of the grain. Never- 
theless, this averaging happens on time scale of the order 
of ergodic time r crg . If we are interested in the properties 
of the spectrum at energy much smaller than the Thou- 
less energy 2?Th (the latter is equal to the inverse ergodic 
time, ETh — T~g), we can assume that the motion of 
electron- hole pairs constituting A^(p, 0), is described by 
equation, averaged over irrelevant coordinates: 



d_ 



2 dz 



0. (83) 



wave packet. As we already discussed, the term d 2 /dz 2 
represents fluctuations of the La ypuno v exponent and is 
important at large energy e ~ ^/A/t e . Below we neglect 
this term. 

The function A^(e, (,z) is defined at < ( < I, where 
I is the length of a trajectory, and has the boundary 
condition A^(e,l,z) = 1 at £ = I, see Table I. The 
other boundary condition exists at z = Inb/L, which 
represents the statistical independence of electron-hole 
wave functions, separated by a distance b. We say about 
these functions that they correspond to a zero mode a. 
Thus the boundary condition is A^(e, (, z = lnb/L) — 
a k , see Fig. 9. 

We write the solution to Eq. (83) as 

A+(e,C,z) = a k e 2ik ^- lnb ^ x 0( j~^ b/L -( + 1 



+ e 2ik ^e{c-i- z ~ x l b,L ) 



(84a) 



Repeating the above analysis for A k (s,x, z), we obtain 



+ e 2 ^e(c- z - ln x b/L ^. 



(84b) 



The density of states is described by a solution of 
Eilcnbcrger equation, represented in the form of series 
in a±(e, x). According to Eq. (58), the quantum diffrac- 
tion leads to a finite separation between points Xj in the 
expression for a±(e, x). The separation is determined by 
/d m in and c/> m i n defined by Eq. (57). We conclude, that 
the Green function may be represented as a sum of func- 
tions A^(s,^,lnp m i n /L) over all positive integers k, see 
Eqs (58) and (77). 

The functions 



a+(s,C) 
a-(s,() 



exp(2«e(C - I)), I - t e < C < h 
ao, < x < I — te, 

exp(-2ie(), < C < t e ; 
ao, t e < C < I- 



(85a) 
(85b) 



Here £ stands for a coordinate of the center of mass of a 



allow us to rewrite the sum for Q(s,Q in terms of 
A^(e, C, z) in a more compact form of Eq. (15). We in- 
troduced a new notation ao for the zero mode, related to 
a by the following equation The functions 

a±(e, x) have all properties we discussed above. They 
satisfy the boundary conditions at one terminal of a clas- 
sical trajectory see Table I, and coincide with semiclas- 
sical solutions of the Eilenberger equation at small dis- 
tance C ^ te or I — ( < TEirom this terminal, when the 
quantum diffraction effects may be disregarded. As an 
electron travels for time interval equal to the Ehrenfcst 
time, te, the quantum scattering suppresses the oscilla- 
tions and a±(e, () jumps to a zero mode value ao- 

Using Eqs. (85) to calculate the density of states at 
large energy, e ^> j csc , we obtain the result of Eq. (82). 
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Indeed, by direct averaging of a±(e, x) in the form given 
by Eqs. (85), and comparing the result with Eq. (62), 
we find a — e~ 2t£TE -f csc /(2ie), and then again integrat- 
ing «4(e, x) over the phase space we reproduce Eq. (80), 
apart of the term A 2 £ 2 te/A 2 , which is important only at 
the extremely large energy e > a/A/t e . 

Below we perform the self-consistent calculations of the 
Green function in the Riccati parameterization, Eq. (15), 



with functions a±(e,x) given by Eqs. (85). Within the 
anzats of Eqs. (85), there exists only one parameter ao, 
which completely determines the Green function in the 
grain. This single parameter can be found with the help 
of Eq. (18). Indeed, we substitute g(e,x) and f(e, x) 
in terms of the Riccati variables a±(e,x), write ao = 
exp(— 2itp) and obtain the following equation for ip: 



tan ip 



- e 



7o SC V cos 2^ 

2r E 



Tesc T E \ COS it 

+ 2(1 - e ~ 7cscTE ) cos 2ip + 2e~ 7racTE sin 2ip In ■ 



cos(f/> — £Te) 



7-cO-nO f sin £ (2r E -Q 



cos ip 



TE 



7csc e — : - + 2sin2^1n- 

L cos el cos{tp + e{l — te)) 



dl 



Representing £(e,x) in terms of functions a±(s, (), we find the density of states from Eq. (8): 

' NT N* 



Si 



7T 2 7 2 
" /esc 



+ e~ 27racTE Im <^ tan 2tp + In 



(n + 1 /2)e- 7r7 -( n+1 / 2 )/ £ + ( 2t W^ + (« + V 2 )) e -^-( n+1 / 2 )/' 

v n=0 n=7V*+l 

r2r E 



COS('0 — £Te) 
C0S(V> + £T E ) 



7e S ce-^Im{ln C °^- £ ;;-^h ^ 
[ cos(V> + e(« - te)) 



(87) 



Here iV* = int(jerE/7r — 1/2) is the integer part of 
jeT E /n — 1/2. The terms in the first line of Eq. (87) 
originate from purely ballistic contributions of the elec- 
tron Green function, while the term, containing tanh 2tp, 
represents the contribution of the zero mode. The re- 
maining terms describe the contribution from the parts 
of the phase space, where both the zero mode and bal- 
listic solutions for the electron-hole wave pairs are inter- 
mixed. In the semiclassical limit only the first term in 
the first line of Eq. (87) survives with N* — > oo, and thus 
coincides with Eq. (52). In the quantum limit, te — > 0, 
only the term containing tanh 2-0 remains in Eq. (87), 
and corresponds to the random matrix result of Eq. (37). 

According to Eq. (87), the density of states is always 
finite at e > 7t/4te, and is formed by at least the second 
term of the first line of Eq. (82). We investigate if the 
gap energy E g is smaller than 7t/4te, so that the density 
of states originates from electron states described by the 
zero mode. 

From Eq. (87) we conclude that the density of states 
also exists below tt/Ate as long as the solution to Eq. (86) 
has an imaginary part. For small energies e, Eq. (86) 
has real solutions and Af(e) — 0. As energy e increases, 
two solutions approach each other and at certain value 
of energy e = E s , they are both equal to some value ip c . 
The energy E g is the gap energy and the density of states 
has a square root singularity just above the gap. Indeed, 
we can expand both sides of Eq. (86) to the lowest order 
in ip — tp c and e — E g . The expansion in e — E g linear, 
but expansion mip — tp c is quadratic, since the first order 



derivative of Eq. (86) with respect to ip vanishes at ip c , 
see Sec. Ill B 1 for a similar analysis. 

Figure 10 demonstrates dependence of the energy gap 
E g on the Ehrenfest time te, shown by a solid line. To the 
lowest order in 7 C scTe the gap energy decreases linearly 
from the random matrix result, 20 see also Eqs. (38) and 
(39): 

E s w ^ mt (l - 0.23 7eS cT E ). (88) 
The density of states above E g is given by 

where c ss c rmt (l — 0.427 csc te). As the Ehrenfest time 
increases, the energy gap decreases and is given by the 
following asymptote at te7 csc 3> 1: 

e e » -T- (l —) ■ (90) 

4r E V T E 7osc/ 

We plot the asymptote of Eq. (90) by a dotted line in 
Fig. 10. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we considered the effect of quantum dis- 
order and quantum chaos on electron density of states in 
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FIG. 10: The gap energy E s as a function of the Ehrenfest 
time te is shown by a solid line. The dotted line represents 
the asymptote result for long Ehrenfest time, te7 CS c 3> 1, 
given by Eq. (90) and the dashed straight corresponds to the 
short Ehrenfest time limit, te7 CS c <C 1, see Eq. (88). 



an Andreev billiard - a metal grain or a semiconductor 
dot of irregular shape, connected to a superconductor. 

First we restate our assumptions about the considered 
system and briefly review prior publications, discussing 
related systems. We study the density of states at en- 
ergies much smaller than the gap energy in the super- 
conductor and the Thouless energy i/ph (or the inverse 
crgodic time r crg ). The opposite limit, when the escape 
time is equal or less than the ergodic time, was consid- 
ered in refs. [34] for Andreev billiard, in refs. [30,31,32] 
for one dimensional systems (wires) and in refs. [5,35] for 
bulk superconductors. 

In our calculations we use the exponential distribution 
P(l) of trajectory lengths I, described by Eq. (47). An 
analysis of motion in classically chaotic systems shows de- 
viation of the distribution function from the exponential 
law due to anomalously long trajectories. 36 The effect 
of these trajectories on the density of states of an An- 
dreev billiard was studied in ref. [37] within scmiclassical 
theory, when the contribution of long trajectories is es- 
pecially important. We assume that the deviation from 
the exponential distribution occurs at lengths I longer 
than either the Ehrenfest time or the mean free time, 
and we neglected the contribution of anomalously long 
trajectories to the electron density of states. 

We study only non-integrable systems. The difference 
between integrable and chaotic non-integrable systems 
was examined for a quantum dot in refs. [20,38]. 

We do not consider the effect of finite mean level spac- 
ing in the grain. The finite mean level spacing leads to 
mesoscopic fluctuations of the energy levels, and particu- 
larly, to the fluctuations of the gap. 39 The fluctuations of 
the energy gap smear the square root singularity of the 
density of states, Eq. (38b), and produce an exponen- 
tially small tail of the ensemble average density of states 
below the gap energy, see refs. [31,39,40,41]. Neverthe- 
less, the corresponding scale for the gap fluctuations and 



the mean level spacing near the gap are determined by 
the following small energy scale: (iJ^osc) 1 ^ 3 -C 7csc- Sig- 
nificant fluctuations of the first energy level (energy gap) 
occur in exponentially rare cases. 31,39 ' 40,41 

We assume that the reflection at the contact with the 
superconductor is described by the Andreev mechanism, 
i.e. the contact has no tunnel barrier between the grain 
and the superconductor. The effect of finite barrier was 
studied in ref. [20] in the random matrix limit, where it 
was shown that a weak (normal) reflection at the barrier 
does not qualitatively change properties of the curve of 
the density of states. 

We do not take into account the Coulomb interaction 
between electrons in the grain responsible for the effects 
of zero bias anomaly and the Coulomb blockade. In our 
system, the conductance of the SN contact is large (the 
number of channels N » 1), and the Coulomb interac- 
tion may be disregarded. The influence of the zero bias 
anomaly on proximity effect was studied in ref. [42] , and 
the Coulomb blockade in a small metal grain connected 
to a superconductor was studied in ref. [43]. 

we also assume that the repulsion in the Cooper chan- 
nel is sufficiently small and we disregard it. The effect of 
repulsion was studied, for example in refs. [44,45,46]. 

The model of clean chaotic Andreev billiards, consid- 
ered in Sec. IV of our paper, was also recently addressed 
in paper [24]. According to ref. [24], the energy gap and 
the density of states above the gap are determined by the 
random matrix theory result 20 with the appropriately 
modified number of channels and the mean level spac- 
ing. Although the results 24 qualitatively coincide with 
ours, they are quantitatively differ from our results, pre- 
sented in Eqs. (88)-(90) and in Fig. 10. We believe that 
the calculations 20 within random matrix theory are jus- 
tified only as te — > 0, since these calculations take into 
account only the zero mode. In general, as we demon- 
strated in Sec. IV, the ballistic parts of trajectories at 
length smaller than te, also known as Hikami boxes, give 
rise to a substantial contribution to both the density of 
states and the self-consistency equation. Consequently, 
Eqs. (86) and (87) cannot be reduced by a change of mean 
level spacing and the number of channels to the random 
matrix theory equations of ref. [20], see also Eqs. (36) 
and (37) in Sec. IIIB 1, and thus the quantitative result 
of ref. [24] is questionable. 

The most interesting qualitative result of our work is 
the oscillations of the density of states as a function of 
£te- [An oscillating character of the density of states 
near the gap was also found in ref. [23].] These oscilla- 
tions are related to the contribution to Af(e) from the 
ballistic part of the electron Green function. Although 
the full length of different trajectories significantly varies, 
the Ehrenfest time is nearly the same for all long trajecto- 
ries with length I > te. Indeed, even strong fluctuations 
of the Lyapunov exponent are averaged out on time scale 
of the ergodic time T erg , resulting in weak fluctuations of 
the Ehrenfest time. The fluctuations of the Ehrenfest 
time are important only at parametrically large energies 
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e oc ^/ A/te and lead to suppression of the oscillations. 
Similar oscillations also exist in frequency dependence of 
the weak localization correction to the conductivity 7 and 
the correlation function of electron energy levels. 8 

We notice that the oscillations discussed above should 
not be associated with the oscillations of the energy lev- 
els in the Wigner-Dyson statistics. The latter oscillations 
have a much smaller period (equal to the mean level 
spacing) than . The energy level oscillations were 
studied in systems with quantum disorder 47 and with 
quantum chaos. 48,49 Similar oscillations are also present 
in the density of states near the gap with the period 
oc <5^ 3 7e/ c 3 . 31 ' 39 In the present paper we do not take into 
account these oscillations, as well as other mesoscopic 
effects. 

The observation of oscillations in the density of states 
with period oc r^ 1 is the major distinguishing charac- 
teristic of quantum chaotic systems from disordered sys- 
tems, since in quantum disorder such oscillations are ab- 
sent. The period of the oscillations may be used to esti- 
mate the Ehrenfest time and the Lyapunov exponent in 
each particular metal grain. 

In conclusion, our results allow us to describe the 
classical-to-quantum crossover in the density of states 
of Andreev billiards. Methods, developed in this paper 
should also be applicable to study this crossover in many 
other problems of solid state physics and optics. 
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APPENDIX A 

In this appendix we show how the small angle scatter- 
ing term, which mimics the effect of quantum diffraction, 
may be used to obtain the results of Sec. IV B 1. 

A small angle scattering potential, with a typical scat- 
tering angle m ; n ~ y/Xp/L, leads to the following form 
of the self energy part of the Eilenberger equation: 

S(e,r,n) = -Q{e, r, n) + ^-W 2 J(e, r, n), (Al) 



TO 



2r Q 



where we introduced the full mean free path r , 
1 



dn 



(A2a) 



and the transport time r q , 

— = 27rzmi / |^(n-n')| 2 |n x n'| 2 ^. (A2b) 

Tq J lid 



The first term in the right hand side of Eq. (Al) does 
not contribute to the scattering term of the Eilenberger 
equation. The second term represents small angle diffrac- 
tion and the Eilenberger equation acquires the form with 
the Perron-Frobenius differential operator. For the two 
dimensional electron system V n = d/d<fi, where <j> is the 
angle in the direction of n. we estimate r q and t below. 

We use the Ricatti parameterization of the Green func- 
tion, see Eq. (15). For the scattering term of the Eilen- 
berger equation in the form of Eq. (Al), we rewrite 
Eqs. (16) as: 



(2ie — C) a + 
(2ie + £)a- 



1 ( r~, . 

1 



V 2 a_- 



2a_ [V n a+] 2 
1 + a + a_ 

2a+ [V n a_] 2 
1 + a + a_ 



(A3a) 
(A3b) 



Here again we omitted variables {e, x} of the functions 
a± = a±(e, x). 

We notice that the average of the electron-hole wave 
function over the phase space is small at large energy 
£ ^> 7csc, see Eq. (62). This observation allows us to 
neglect non- linear terms in Eqs. (A3), since these terms 
would contribute to the higher than the second order in 
7c SC /e. We write down equation for A^(e, xi, X2) using 
the linearized equations of motion for a±(e, x), Eqs. (A3). 

Ate T£ cm -£(p, 4>) + ^l Af(e, R, N, p, 4>) = 0,(A4) 
2r q 

where £ cm and C(p,4>) are defined by Eqs. (65). As we 
already mentioned, the non-linear terms in the right hand 
side of Eqs. (A3) lead to higher terms in 7 CS c/e- 

After averaging over the position of the center of mass 
in the phase space, we obtain a counterpart of Eq. (66) : 



Aie + s 2 7csc - C, 



V 2 , U 2 t (e,p,(A)=7esc. (A5) 



2r q 

Below we assume that s 2 ~ 1. Introducing new variables 
z and x, according to Eq. (67), and then averaging over 
X, we arrive 7 to the following equation for Af(uj, z) with 
u> = Ae - ij csc : 

ilu + X8 Z + ^dl + ^ ( -i) 1 + 7 9 z 



2r a 



^2 Tesc 

(A6) 

where d z — d/dz, and 7 is a numerical coefficient of 
order of unity. 7 Equation (A6) should be supplied with 
the boundary conditions at z = \nb/L, determined by 
Eq. (69). The solution can be found as a superposition 
of two terms, see Eq. (70). One term, A^ b (e, z), is given 
by Eq. (71) and the other term is given by 



A ± 



2.1: 



00 



w(u>, z) 



;-(A7) 



w(lj, In b/L) 

Function w(u, z) is a solution 7 of homogeneous Eq. (A6): 



w(ui, z) = exp 



/ iuj 


X 2 uj 2 \ 


[Tx- 


AX 3 ) 



In 



Xr n 



Xr q e 2z + 7/2 



(A8) 
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Now solutions for A 2 (e, z) are regular functions as z — > 
— oo, and correspondingly the product of two electron- 
hole pairs, see Eq. (63), can be found at coincident points. 

An exact form and strength of small angle scattering 
potential is not important, since as we will show later, 
the density of states depends on the strength of the po- 
tential only as a logarithm of the ratio of the transport 
time of this potential and the size of the SN contact. 
We will choose the strength of the scattering potential 
so that at small time this potential scattering would de- 
scribe the quantum mechanical diffraction effect. The 
scattering time r q may be estimated as r q = To^in' 
compare Eqs. ( A2a) and (A2b) . We assume that Ato ~ 1 
(since both A and To are characterized by the same poten- 
tial energy of electrons in the dot) and then we estimate 
Ar q ~ 0~ 2 n = L/X F , see Eqs. (A2a) and (A2b). Combin- 
ing Eqs. (69) and (A7) and we obtain Eq. (76). 

In conclusion, we discuss the meaning of Eq. (A8). 
Within logarithmic accuracy, the term Ar q e 2z +7/2 can 
be replaced by max{Ar q e 2z , 7/2}. At large scales p and 
of a wave packet, z > z c = In ^/2Ar q /7, we can disre- 
gard the quantum diffraction, described by the last term 
in Eq. (A4). In this "classical" approximation we find 







)ln 2 ^] 




4A 3 , 





At smaller scales of the wave packet, z < z c or p 2 + 
L 2 (f> 2 < Pmin' the quantum diffraction overpowers the 
Lyapunov classical divergence of trajectories. Conse- 
quently, the function w(e, z) weakly depends on the scale 
z, see Eq. (A8), Eq. (A8), and we use w(e, z — > —00) 
w(e, z c ). It is exactly this property of the function w(e, z) 
was used in Sec. IV Bl of the present paper, when we 
used w(e, z c ) 



In this Appendix we used linearized Eqs. (A3), which 
are justified at high energy e ^ 7 CS c- In general, func- 
tions Af(e,z) are described by non-linear Eqs. (A3). 
Nevertheless, the right hand side of Eqs. (A3), which 
contains non-linear terms, is necessary only to determine 
the function w(e, z) at small scale of the wave packet, i.e. 
for z < z c , and at already at scale z ~ z c the right hand 
side of Eqs. (A3) does not affect the solution for w(e, z) 
within logarithmic accuracy. That is why the non-linear 
term is not important within logarithmic accuracy. 



Within the same logarithmic accuracy, we use approx- 
imate solutions for A^(e,z) at z ~ z c for k > 2 in 
Sec. IV B. 
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